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We introduce a new mesoscopic model for nematic liquid crystals (LCs). We extend the particle- 
based stochastic rotation dynamics method, which reproduces the Navier-Stokes equation, to 
anisotropic fluids by including a simplified Ericksen-Leslie formulation of nematodynamics. We 
verify the applicability of this hybrid model by studying the equilibrium isotropic-nematic phase 
transition and nonequilibrium problems, such as the dynamics of topological defects, and the rhe¬ 
ology of sheared LCs. Our simulation results show that this hybrid model captures many essential 
aspects of LC physics at the mesoscopic scale, while preserving microscopic thermal fluctuations. 
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I. INTRODUCTION 

Liquid crystals (LCs) possess anisotropic interactions 
because of their molecular shape. For example, molecules 
with a rod-like rigid core tend to align parallel to each 
other and form a mesophase called nematic; molecules 
with a disk-like rigid core form discotic phases. In both 
cases the rigidity is generated by different combinations 
of aromatic rings [1]. Macroscopically, this anisotropy 
leads to a series of phase transitions that break rotational 
and translational symmetries in a step-wise fashion. Be¬ 
cause of their capacity to reorient, also in response to 
external fields, LCs are used in a wide range of applica¬ 
tions: from the ubiquitous electronic displays, to micro¬ 
lasers [2-4] and lubricants [5]; but they are also rising 
to an important role in biomedical sciences and applica¬ 
tions [6], and in our comprehension of morphogenesis and 
evolution of living organisms [7]. 

Hydrodynamic flow can also couple with the local 
preferential direction (director) established in nematic 
LCs [8]. Recently, Sengupta et al. have extended mi¬ 
crofluidic applications to anisotropic fluids and have 
found surprising topologies of the orientational field 
[9, 10], and explored fluid and colloidal transport [II- 
13]. These effects originate from the intimate connec¬ 
tion between LC rheological properties [14] and their lo¬ 
cal alignment, which can be easily controlled. However, 
the complex interplay of confinement to a mesoscopic 
scale, i.e. of the order of /rm, surface interactions, hy¬ 
drodynamic flow and generation of topological defects 
still poses formidable challenges to both theoretical and 
experimental investigations. 

The theoretical description of LCs based on static con¬ 
tinuum theory started in the early 1920s with the work 
of Oseen [15], Zocher [16], and Frank [17]. The earliest 
dynamics theory of LCs can be dated back to 1931 by 
Anzelius [18]. In the 1960s Ericksen and Leslie devel¬ 
oped a hydrodynamics theory [19-23] based on the LG’s 
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velocity field v{r) and a unit vector describing the lo¬ 
cal director d{r). The nematodynamic equations of the 
Ericksen-Leslie model are widely used but rest on the as¬ 
sumption that the nematic order parameter is a constant 
and the nematic LC is uniaxial, and therefore they can¬ 
not describe physical situations where there is a strong 
variation of the the nematic order, such as the isotropic- 
nematic phase transition and the dynamic of topologi¬ 
cal defects. For cases where there is a strong variation 
of the nematic order a tensorial description is necessary, 
such as the Beris-Edwards formulation [24] or the Qian- 
Sheng formulation [25]. These two approaches differ in 
the form of the elastic free energy, that is, the former con¬ 
siders the elastic free energy in the one-constant approx¬ 
imation, while the latter allows for two different elastic 
constants. We are not aware of any tensorial description 
of nematodynamics in terms of all three elastic constants, 
i.e. splay, twist and bend. 

Microscopic models have in general the advantage of 
providing detailed dynamics at small spatial and tem¬ 
poral scales. Molecular dynamics simulations are well 
suited for this task. However, physical phenomena at 
the mesoscopic scale are still so computational demand¬ 
ing that they are out of the reach of atomistic simula¬ 
tions. This predicament can be ameliorated by adopting 
a coarse-grained description of the molecular degrees of 
freedom that effectively includes hydrodynamic modes. 
One widely-used method is the lattice Boltzmann scheme 
that simulates the evolution of the Boltzmann equation 
for a simple fluid on a regular lattice by a series of colli¬ 
sion and propagation steps for probability density func¬ 
tions defined on the sites of a lattice [26]. This scheme 
has been generalized to LCs with successful results [27- 
29]. However, lattice Boltzmann schemes have the lim¬ 
itation that they do not include thermal fluctuations. 
A different approach that emphasize the particle aspect 
is the dissipative particle dynamics (DPD). This is an 
off-lattice, particle-based method, where particles repre¬ 
sents fluid elements and are subject to pair-wise additive 
forces, which conserve momentum locally and thus gen¬ 
erate the Navier-Stokes hydrodynamics for simple fluids. 
The DPD scheme has also been extended to LCs [30, 31]. 
Whereas this method has been shown to reproduce equi- 



2 


librium phase diagrams of LCs [30, 31], we are not aware 
of any attempt at reproducing nematodynamic behavior, 
though the original version of DPD does reproduce the 
correct hydrodynamics of simple fluids [32], 

More recently, Malevanets and Kapral [33] introduced 
the stochastic rotation dynamics (SRD) model. This is 
also an off-lattice, particle-based model where each par¬ 
ticle represents a small parcel of fluid. The fluid evolves 
through a series of collisions and streaming steps that ex¬ 
actly conserve mass, linear momentum and energy, and 
additionally respect Galilean symmetry. Thus, the cor¬ 
rect hydrodynamic modes are generated. Because SRD 
is a particle-based method, fluctuations are naturally 
present. The SRD method has been applied to a variety 
of systems, from colloids [34] to polymers [35, 36], from 
the modeling of the solvent carrying hydrodynamic inter¬ 
actions in vescicle self-assembly [37] to the flow-induced 
shape transition in red blood cells [38]. In general the 
SRD model is useful whenever both thermal fluctuations 
and hydrodynamic modes are physically important. 

In the present work we extend the SRD model to 
anisotropic fluids, namely LCs. This is primarily done 
by giving orientational degrees of freedom to the SRD 
particles which obey the (simplified) equations of nema- 
todynamics in the formulation of Ericksen-Leslie. For the 
sake of simplicity we restrict ourselves to methods and 
results in two dimensions (2D). We show below that the 
SRD is amenable of studying both the equilibrium and 
nonequilibrium behavior of LCs. This new extension of 
the SRD model opens the door to the investigation of a 
wealth of LC phenomena occurring at the mesoscale. 

This work is organized as follows. In Sec. II we describe 
the theoretical background of the equations used in the 
present model, and the details of the numerical imple¬ 
mentation. We validate this model for nematic LCs by 
considering three study cases in Sec. Ill; i) the isotropic- 
nematic phase transition; ii) the production and annihi¬ 
lation of topological charges, and iii) LC rheology under 
shear flow. Finally, we summarize and discuss our results 
in Sec. IV. 

II. THE MODEL 
A. Theoretical Background 

We start by considering the standard SRD model that 
has so far been used to describe fluids with isotropic in¬ 
teractions. The system is composed of N particles of 
mass TTii with positions r)(t) and velocities Vi{t), where 
i G [l,.^]- The evolution in time t proceeds through a 
series of two steps: 

(i) the free-streaming step 

ri{t + St) = ri{t) + Vi{t)dt , (I) 

where all positions are updated; 

(ii) the collision step 

Vi(t + St) = uci(t) +'R[viit) - uci (t)] (2) 


where all velocities are updated by rotating the fluctuat¬ 
ing part of the velocity with respect to the center of mass 
velocity 

UC. {t) = = X] 

2=1 i—1 

In Eq. (2-3) the calculations are performed in a cell-wise 
fashion, that is, the system is divided with a regular grid, 
and uCi is computed from the JVci particles within the 
cell Ci to which particle i belongs. 

The rotation (or collision) matrix R is orthogonal, 
R~^ = R^, where the superscript T denotes transpo¬ 
sition. It rotates, independently in each cell, the fluctu¬ 
ating part of the velocity by an angle a about an arbi¬ 
trary axis. Because of its action on the fluctuating part 
of the velocity the collision rule conserves momentum; 
because of its orthogonality the energy is also conserved. 
From these facts follows that the fluid obeys the Navier- 
Stokes equations. Malevanets and Kapral showed [33] 
that if R satisfies detailed balance then the fluid ap¬ 
proaches a Maxwell-Boltzmann distribution and obeys 
the iL-theorem. In practice there are many ways to 
choose R so that all the requirements are met. In 2D 
R can only rotate the velocities by an angle ±a with 
equal probabilities; in 3D it is common, e.g., to choose 
a random axis and rotate around it of a fixed angle. To 
ensure Galilean invariance the grid must be shifted by 
a random amount at each time-step [39] so that artifi¬ 
cial correlations among particles do not build up due to 
repeated collisions with the same neighbors. 

We now come to the extension to nematic LCs. The 
particle’s degrees of freedom must be augmented by a 
unit vector describing its orientation di, = 1. Two 

particles separated by a distance ||r) — r)|| ^ e interact 
through the Lebwohl-Lasher potential [40] 

D. = - ^(d* • 4)2. 

(iJ) 

where (i, j) indicates that i and j are neighbors. 

Lin et al. [41, 42] proposed a simplified version of the 
Ericksen-Leslie equations. These are nonparabolic, dis¬ 
sipative equations that describe the flow of nematics in 
the one-constant approximation. Although they repre¬ 
sent a drastic simplification of the original Ericksen-Leslie 
model, they retain the essential mathematical features, 
and they obey an energy law similar to the one used in 
the Ericksen-Leslie model [42]. 

We propose the following hybrid model of SRD for 
anisotropic fluids 

dv 

— + V ■ Vv = V ■ (vVv) — VP/p — W ■ TV (4) 

f)(] —> —> —+—*—> 

_ + . Vd - d . Vu = 7.. V^d - 7/(d) + m (5) 

where p is the density, P the pressure, v = r}/p the kine¬ 
matic viscosity, and tt = (Vd ■ Vd) the Ericksen-Leslie 
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stress tensor, that is the tensor whose (q;/ 3) component is 
dd/dva ■ dd/drp. In Eq. (5), 7 ^,^ is the elastic relaxation 
constant, the term f{di) = dUi/ddi is the molecular field 
inducing nematic ordering, 7 its strength, and ^(t) is a 
Gaussian white noise for the director’s angular velocity, 
with {^a) = 0 and = 2k^T-i5afi5{t - t'), 

where fce is Boltzmann’s constant, T is the temperature, 
5ap is a Kronecker delta, and 5{t) is Dirac’s delta distri¬ 
bution. 

The Ericksen-Leslie stress tensor tt represents the feed¬ 
back of the director field d{r, t) to the bulk flow of 
molecules. In this formulation, tt is directly responsible 
for the non-Newtonian behavior of the LC flow. Because 
of their molecular anisotropy the viscosity of a LC fluid 
does not depend solely on the shear stress r = rjdv/dx, 
but the director field d{f, t) also participates in the vis¬ 
cosity generation. This means that the velocity held can 
couple to the director and reorient it, and, also, that a 
reorientation of the director may generate a how, usually 
called backflow. 


B. Numerical Implementation 


We now describe the numerical implementation of our 
hybrid model for LCs. We will make a distinction be¬ 
tween quantities computed on a particle level, such as 
Vi, and on a cell level, such as uCf We note that the 
standard SRD steps described in Eqs. (1-2) recover the 
Navier-Stokes equation, that is Eq. (4) without the last 
term. Also, Eq. (5) applies in the Eulerian picture (lab- 
frame) to the huid parcel moving along the streamline. 
Since in simulations the director of the ith particle, di, 
moves along the how, the convective term v ■ Vd is ab¬ 
sorbed in this Lagrangian picture (comoving frame). 

The SRD algorithm for LCs consists of the following 
steps: 

(i) free streaming 

n{t + 5t) = ri{t) + Vi{t)5t, (6) 

this is identical to the standard SRD step; 

(ii) cell-wise calculations, that is, the particles are 
grouped in different cells and uci, dci, gradients such 
as Vd and Vu, and the Ericksen-Leslie elasticity tensor 
TT are also calculated; 

(hi) LC alignment 


di(t -f 5t) — di(t) -\- 


d-Vv + -f^^V^d 


St (7) 


where Eq. (5) is implemented. 

(iv) Collisions and backhow. The cell-wise, center of mass 
velocity is calculated as in Eq. (3), and the contribution 


from the Ericksen-Leslie tensor is then added 

dci (^) = {t) + AV • TTCi (8) 

v^(t + St) = uc^{t) + AhR. [vi{t) - uci{t)] (9) 

We consider here only a 2D system, thus the matrix R 
rotates the particles’ thermal velocities around one axis, 
conventionally denoted as the z-axis, either clockwise or 
counter-clockwise by a hxed angle a stochastically. The 
parameter Ah is the thermostat scaling factor whose role 
will be explained below. 

In general, SRD does not conserve angular momen¬ 
tum. To impose the conservation of the fluid’s angular 
momentum in 2D the angle a must be chosen [43, 44] 
such that 

2AB 

= - A^ + B^ ’ A2+R2 (10) 


where 


J^Ci 

^ = {v^-ucj\\z, (11a) 

^fci 

B = '^fi ■ {vi - uci) ■ (11b) 

i=l 

A thermostat for translational and rotational velocity 
is implemented to control the temperature of the system. 
Equipartition of the energy is applied to keep the same 
amount of energy in each degree of freedom. The temper¬ 
ature is defined as k^T = ^ ^ - ^rriivf from the trans¬ 
lational kinetic energy, or as feeT = ^ from the 

rotational kinetic energy in 2D, where Ii is the moment 
of inertia of the particles and Wi is the angular velocity. 
We employ a simple velocity rescaling thermostat that 
scales linear and a ngular v elocities in a cell-wise fashion 
by a factor Ah = i/T/Tci (see Eq. (9)), where is the 
instantaneous kinetic temperature in cell Ci. 

We use two types of boundary conditions (BC), peri¬ 
odic BC in both the x and y directions for simulations 
of bulk systems, and a no-slip wall perpendicular to the 
X axis (while periodic BC are implemented in the other 
direction) for simulations of a channel geometry. The 
no-slip BC are implemented with the usual bounce-back 
rule, that is, the velocity of the LC particle is inverted 
upon collision with the solid wall 

^new ^ _^old (^2) 

Additionally, “ghost” particles [45] are used at the walls 
to fill partially occupied cells up to the average particle 
density. At solid interfaces the anchoring, that is, the 
preferential angle between LC particles and walls needs to 
be specified. Homeotropic anchoring is implemented by 
placing ghost particles in the walls with their orientations 
aligned perpendicularly to the walls. 

The system is initialized by assigning the positions 77 
and velocities A- The particles are uniformly distributed 







in space. The linear and angular velocities are assigned so 
that they are distributed according to a Maxwell distri¬ 
bution at the target temperature, and that the equipar- 
tition theorem is obeyed. 

We perform all our simulations with the following pa¬ 
rameters: the SRD rotation angle is a = 120°, particle 
mass rrii = 1, moment of inertia /^ = 1, and the elas¬ 
tic and molecular relaxation constants are 7 ^,^ = 10“^ 
and 7 = 8 X 10“"^. It is a common choice to set the 
SRD grid size Sx = 1 and the timestep St = 1. At 
each time step the grid is shifted by a random displace¬ 
ment vector with components uniformly distributed in 
the interval \—5x/2,5x/2]. In the following we mea¬ 
sure the temperature T in units of where 

Vmax = SxjSt is the maximum propagation speed of a 
particle (related to the Courant condition). The mean 
free path A = St^k^aTjmi is typically smaller than one 
but the grid-shift method avoids the build-up of spurious 
correlations. 


III. APPLICATIONS 
A. Isotropic-Nematic Phase Transition 

We perform simulations of a system in equilibrium at 
fixed T with periodic BC in both x and y directions. 
A useful way to characterize the degree of ordering of a 
nematic LC is to define an order parameter that is zero 
in the isotropic phase and one in the nematic phase. It is 
common to take the largest eigenvalue S of the nematic 
order tensor 

Q = (^2^* 27 ’ 

which is a traceless, second-order tensor, and where ® is 
the dyadic product and I is the unit tensor. 

We consider a system of iV = 150 000 particles, and 
of size Lx = Ly = 50, subdivided with a 50 x 50 grid, 
thus the average number of particles per cell (AcJ = 60. 
Figure 1 shows the dependence of the nematic order pa¬ 
rameter S on T, at fixed density. At low T the ne¬ 
matic order parameter is very close to one, indicating 
nearly perfect nematic alignment. As T increases S de¬ 
creases continuously until it reaches the constant value 
1/4, which characterizes the isotropic phase in 2D [46]. 
We find that the isotropic-nematic transition occurs at 
T ~ 0.18. A LC system in two dimensions or higher 
undergoes a temperature-driven isotropic-nematic phase 
transition. This phase transition is continuous in 2D, 
as expected for a system with the symmetry of the XY 
model (though there is some controversy [47-49]). 

Contrary to the hydrodynamic Ericksen-Leslie model, 
our SRD hybrid model can reproduce the isotropic- 
nematic phase transition. The reason for that is the par¬ 
ticle nature of our model. Different parts of the system 



FIG. 1. Nematic order parameter as a function of thermostat 
temperature. The solid line is just a guide for the eye. A 
continuous isotropic-nematic phase transition occurs at T ~ 
0.18. 

may have different local orientations and thus, at high 
enough T, the system is globally disordered. 


B. Dynamics of Topological Defects 

During a quench from the isotropic to the nematic 
phase a LC substance develops many topological de¬ 
fects [50] which influence the kinetics of the phase tran¬ 
sition. Solid boundaries imposing geometric or ener¬ 
getic (anchoring) constraints, or external fields in gen¬ 
eral may induce different, local orientations in a nematic 
LC. Because of the conflict between these local orienta¬ 
tions topological defects are produced. These are regions 
of the fluid where the local nematic order parameter is 
vanishing and the director field is undefined [51]. 

Although there were attempts of using Ericksen-Leslie 
model to simulate topological charges [52], those struc¬ 
tures cannot be properly treated since the nematic order 
parameter is essentially constant in that approach. How¬ 
ever in our model, the nematic order parameter in the 
macroscopic scale is calculated from the mesoscopic SRD 
particle configurations, hence the formation of topolog¬ 
ical defects is an essential feature in this particle-based 
approach. 

We investigate the generation of topological defects 
and their dynamics by quenching the system from a 
high value of T, well in the isotropic phase, into the ne¬ 
matic phase at T = 0.053, with an equilibrium value of 
S = 0.91. Figure 2 shows the temporal evolution of the 
topological defects formed during the quench. A large 
number of ±1/2 charges are generated. The total topo¬ 
logical charge of the system is zero, and, by conserva¬ 
tion of charge, an equal number of positive and negative 
charges are present. Immediately after the initial quench 
the SRD particles start to align with their neighbors, due 
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FIG. 2. Snapshots of the system at different times showing the 
temporal evolution of local nematic order parameter (color) 
and local director (black dashes). As the system is quenched 
from the isotropic to the nematic phase topological defects 
are formed. In panel (d) the red dashed circle marks the q = 
— 1/2 topological defect, while the red solid circle indicates 
the q = +1/2 defect. 


to the torque induced by the molecular field /(d). After 
some time the particles start to form larger nematic do¬ 
mains. At this stage (Fig. 2(a)), there are many topo¬ 
logical defects, carrying topological charges of g = +1/2 
and q = —1/2. At a later stage, those charges with op¬ 
posite signs are attracted to each others and eventually 
annihilated with each others (Fig. 2(b-d)), leaving the 
area charge-free as the initial configuration. 

The speeds of topological defects depend on the 
charges they carry, that is, defects with charge q = +1/2 
move faster than q = —1/2 charges; this has been verified 
both by numerical [53-55] and experimental work [56]. 
We find a similar behavior in our simulations. Figure 3(a) 
shows the time dependence of the positions of two topo¬ 
logical charges. Initially, the charges approach each oth¬ 
ers with rather slow velocities, but when the two charges 
are very close they accelerate until they annihilate. It is 
clear the velocity of the positive charge is higher than the 
negative one. The reason for this asymmetry is that the 
velocity of the defect core depends on the sign of the spa¬ 
tial derivative of the director field Vd. Similar result has 
been reported in [54] where lattice Boltzmann simula¬ 
tions of a tensorial formulation of nematodynamics were 
performed. 

To further verify the validity of our model we test the 
time dependence of the separation D between opposite 
charges. Denniston predicted [57] a simple scaling law 




FIG. 3. The position of topological charges as function of 
time is shown in the upper panel. The blue dots (empty cir¬ 
cles) represent the position of positive (negative) charge. The 
red solid (dashed) line are the fits of the positive (negative) 
charge. The lower panel shows their separation as a function 
of time. The blue empty circles are the direct subtractions of 
positions and the red line is the fitting curve using Eq. (14). 


D{t) = c{ta-tY^'^ (14) 

where ta is the annihilation time and c is a constant. 
Figure 3(b) shows that the separation between oppo¬ 
site charges does indeed follow the behavior predicted 
by Eq. (14). Although the scaling law (14) was de¬ 
rived in conditions of no backflow [57], it is still valid 
when the coupling of flow and director field is present, as 
observed in the tensorial treatments of nematodynam¬ 
ics [27]. Then we conclude that our particle-based ap¬ 
proach to nematodynamics can correctly describe regions 
of the fluid with strong gradients of the nematic order pa¬ 
rameter, such as topological defects. 


C. Shear flow 

Shear flow experiments are the canonical way to study 
the rheological properties of fluids. The coupling of the 
elastic deformations of LCs with the transport and de¬ 
formation of fluid elements gives rise to a more com¬ 
plex situation than in the flow of an isotropic fluid. In 
the original formulation of Ericksen and Leslie [19-23] 
the viscous response is characterized by six coefficients 
Oi, i = 1...6, called Leslie coefficients. Later, Parodi 
showed [14] that from Onsager’s reciprocity theorem fol¬ 
lows 06 — ck 5 = ck 2 + 03 , thus the number of independent 
viscosity coefficients is five. 

Our goal in this section is to verify if our hybrid model 
is capable of reproducing the known LC rheology. To 
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FIG. 4. Top row: snapshots of the nematic order parameter (color) and director fields (black dashes) under shear flow are 
shown. The walls at a; = 0 and x — move in the y direction with velocities are Vy = 0.15 and Vy = —0.15 respectively. From 
left to right, panels (a)-(c) show the LC for Ericksen-Leslie stress constant A = 0, 0.05, 0.1, respectively. Bottom row: from 
left to right, panels (d)-(f) show the y-averaged flow velocity as a function of x, Vy{x) for the same values of A as in the top 
row. Shear banding caused by non-zero A is clearly seen by the kinks in flow profiles. 


generate a shear flow in a 2D simulation domain, we con¬ 
sider two no-slip walls at x = 0 and x = and periodic 
BC in the y direction. Because of the presence of solid 
walls the boundary conditions for the director held need 
to be specihed. We employ homeotropic anchoring, that 
is, the LC particles prefer to orient perpendicularly to 
the walls. This is achieved by assigning a perpendicular 
orientation of the ghost particles at the walls. 

The two walls move with equal and opposite speeds 
in the y direction. After some time a stationary shear 
how is generated. A dimensionless measure of shear 
is given by the Weissenberg number W = yr where 
7 = dvyjdx is the shear rate and r is a relaxation time. 
We calculate t from the orientational correlation func¬ 
tion Cl (At) = {d{t + /Si) ■d{t)), where the angle brackets 
indicate ensemble average. We present results for a hxed 
W = 2.04 at T = 0.28 (we note that the presence of walls 
inducing homeotropic alignment shifts the phase diagram 
with respect to Fig. 1, so that T = 0.28 corresponds to 
the nematic state). The system of iV = 75000 particles 
and size = 15, Ly = 50 is divided with a grid of 
15 X 50, thus with {Afci) = 100. 

The steady-state director helds and bulk how prohles 
of sheared LC are shown in Fig. 4 for different Ericksen- 
Leslie stress constant A. For the case of A = 0 the how 
prohle is a straight line as for a Newtonian huid. For 


higher values of A we observe a shear banding effect. By 
increasing A the how prohle starts to have a kink near 
the walls. Shear banding is a nonequilbrium transition 
to a state where regions with diherent shear rates coex¬ 
ist, thus visible through diherent slopes of the velocity 
prohle Vy{x) in Fig. 4. Nematic LCs can produce shear 
bands [58, 59]. The existence of this phenomenon is a 
consequence of the competition between elastic energy 
and viscous dissipation, which in turn is caused by the 
feedback from LC orientation to the how induced by tt. 

The local director held d{r) (represented as black 
dashes in Fig. 4(a)-(c))shows that while the director is 
aligned homeotropically at the walls, thus satisfying the 
BC, it is tilted in the central region of the channel. This 
tilt angle is due to a well-known how alignment mecha¬ 
nism in LCs. The local nematic order map (color code in 
Fig. 4(a)-(c)) shows a subtle change from a high degree of 
nematic alignment at the walls, induced by the anchoring 
conditions, to a slightly lower value in the central region 
of the channel. 


IV. CONCLUSIONS 

We have introduced a new mesoscopic LC model. Our 
model is based on the stochastic rotational dynamics 
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scheme, which is a particle-based algorithm that ignores 
the computationally heavy molecular interactions but 
correctly generates the hydrodynamic modes described 
by the Navier-Stokes equations. 

The model introduced here fills an important gap for 
mesoscopic simulation techniques of LCs. Thermal fluc¬ 
tuations are explicitly present, which are of fundamental 
importance in both thermodynamic and dynamic pro¬ 
cesses, but are neglected in other popular schemes, such 
as the lattice Boltzmann approach. Furthermore, the 
mesoscopic scale and the hydrodynamic behavior is di¬ 
rectly addressed. 

We have shown that this model can be used to study 
various aspects of the physics of liquid crystals. We have 
considered three study cases. First, we have found that 
our model system undergoes an equilibrium phase tran¬ 
sition from nematic to isotropic as the temperature in¬ 
creases. The transition is found to be continuous, as it 
is expected to be in 2D. Second, we have studied the 


nonequilibrium dynamics of topological defects emerging 
in a quenched LC. Topological defects with ±i charge are 
formed but quickly annihilate with each other. The tem¬ 
poral dependence of the distance between two opposite 
charges is found to match remarkably well the theoreti¬ 
cally predicted power law. Third, we have considered a 
shear flow situation and found that the LC system devel¬ 
ops shear bands as the coupling parameter between flow 
and director reorientation is increased. We conclude that 
the model captures the non-Newtonian character of LC 
rheology. 

This hybrid algorithm can be easily applied to complex 
geometries of the confining walls. For the sake of simplic¬ 
ity we restricted the present work to 2D. However, any 
realistic implementation requires a 3D setup. A gener¬ 
alization of the model to 3D is under way and will be 
discussed elsewhere. 

We gratefully acknowledge helpful conversations with 
Martin Brinkmann, Stephan Herminghaus and Thomas 
Hiller. 
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